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Abstract 

A method for performing progressive damage modeling in composite materials and structures based on 
continuum level interfacial displacement discontinuities is presented. The proposed method enables the exponential 
evolution of the interfacial compliance, resulting in unloading of the tractions at the interface after delamination or 
failure occurs. In this paper, the proposed continuum displacement discontinuity model has been used to simulate 
failure within both isotropic and orthotropic materials efficiently and to explore the possibility of predicting the 
crack path, therein. Simulation results obtained from Mode-I and Mode-II fracture compare the proposed approach 
with the cohesive element approach and Virtual Crack Closure Techniques (VCCT) available within the ABAQUS 
(ABAQUS, Inc.) finite element software. Furthermore, an eccentrically loaded 3 -point bend test has been simulated 
with the displacement discontinuity model, and the resulting crack path prediction has been compared with a 
prediction based on the extended finite element model (XFEM) approach. 

I. Introduction 

Carbon fiber-reinforced polymers (CFRPs) are used widely for many applications including aircraft, 
automotive, and civil structures. Some of the appealing properties of CFRPs include high specific strength and 
stiffness, corrosion resistance, fatigue resistance, and suitability for stealth applications. Increased usage of CFRPs 
in the future depends on the development of damage-tolerant composite structures while lowering the cost of 
manufacturing. One of the issues with CFRPs for structural applications is the joining of sections made with fiber- 
reinforced composites. Among other joining methods in current practice, two of the most popular methods are 
adhesive bonding and mechanical fastening where the mechanical fastening is preferred for thicker sections while 
adhesive bonding is used for thinner sections. Delamination is a major type of failure associated with adhesively- 
bonded composite joints. 

In the past, failure of composites has been investigated extensively from the micromechanical and 
macromechanical points of view. On the micromechanical scale, failure mechanisms are related to the properties of 
the constituent phases, i.e., matrix, reinforcement, and the interface. This paper will investigate the advantages and 
disadvantages of modeling damage propagation in materials utilizing a displacement discontinuity approach. The 
displacement discontinuity approach known as ECI (evolving compliant interface) (Refs. 1 and 2) enables 
exponential evolution of the interfacial compliance, resulting in unloading of the tractions at the interface after 
failure/delaminations occurs. The current study utilizes the implementation of this approach available within 
NASA’s FEAMAC (Ref. 3) multiscale micromechanics analysis framework. Given the fact that the commercially 
available techniques such as VCCT (Ref. 4) and cohesive element (Ref. 5) approaches require a priori knowledge of 
the crack path, and as discussed elsewhere (Ref. 6), the current implementation of these approaches present 
difficulties with numerical convergence, forcing one to deal with a number of parameters to overcome these 
convergence issues. Therefore, in this paper, a computationally-efficient multiscale modeling framework is 
presented, and the propagation of damage is investigated utilizing this framework in the context of finite element 
method. The future intent is to apply this approach to the failure/delaminations of composite structures. 

To begin the investigation, VCCT and cohesive element approaches are used to model self-similar cracks in 
double cantilever beam (DCB) the end-notched flexure (ENF) tests, and the results are compared with FEAMAC 
predictions. Once the material data and the simulation-specific parameters have been identified through these 
preliminary simulations, the focus shifts to an eccentrically-loaded three-point bend test. In the absence of 
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experimental data to validate the FEAMAC prediction for the crack path, the extended finite element method 
(XFEM) (Ref. 7) available in the Abaqus finite element package is utilized for comparison purposes. The study 
concludes with highlighting the advantages and disadvantages of the FEAMAC and XFEM approaches. 

II. Displacement Discontinuity Approach and the Evolving Compliance Interface Model 

The displacement discontinuity approach employed herein follows the basic form proposed by Jones and 
Whittier (Ref. 8), in which a debonded interface is modeled by imposing a discontinuity (i.e., “jump”) in a given 
displacement component across the interface. This basic form can be written, 

[ujJ=Rjajf ( 1 ) 


Where [ujf is the discontinuity in displacement component j at interface /, g 7 is the corresponding stress component 
at interface /, and Rj is a proportionality constant that may be thought of as a flexibility for the interface (as it relates 
the interfacial displacement jump to the interfacial stress). This flexible interface model has been applied extensively 
to model fiber-matrix interfacial debonding in composite materials. It was incorporated within the method of cells 
micromechanics model by Aboudi (Ref. 9) and within the generalized method of cell by Sankurathri et al. (Ref. 10). 
Achenbach and Zhu (Ref. 11) added the condition to this flexible interface model that, when in compression, the 
interface must behave as if well-bonded (i.e., Rj = 0). Wilt and Arnold (Ref. 12) added finite bond strength to the 


interface, 



, such that the interface behaves as perfectly bonded until the interfacial stress reaches this 


interfacial strength. Bednarcyk and Arnold (Refs. 1 and 2) further extended the Jones and Whittier (Ref. 8) flexible 
interface model to admit time dependence in the Rj parameter. By allowing this parameter to evolve with time, it is 
possible to enable unloading and redistribution of the interfacial stress upon initiation of debonding or delamination. 
Considering the components normal (n) and tangential ( t ) to a particular interface separately, the equations for the 
displacement discontinuity components can be written as, 


R n (t)o n , 


1 > n 
- g DB 


( 2 ) 


[ u t^ = R t{t)^ t \ I ’ 


- g db 


( 3 ) 


where the time-dependence of the Rj parameter, as well as the finite bond strength of the interface, has been denoted. 
This form of the Jones and Whittier (Ref. 8) flexible interface model has been called the evolving compliant 
interface (ECI) model in order to distinguish it from the previous form that incorporates a constant value for the Rj 
parameter. 

Many forms for the time-dependence of the Rj parameter, both implicit and explicit, were examined by 
Bednarcyk and Arnold (Ref. 1). It was determined that explicit exponential time-dependence provides the model 
with the desired ability to unload interfacial stress, permitting redistribution of this stress. Further, the ECI model 
compares well (in a qualitative sense) with several previous interfacial constitutive models that incorporate 
interfacial stress unloading, e.g., Needleman (Ref. 13). The ECI model representation of the debonding parameter 
put forth by Bednarcyk and Arnold (Refs. 1 and 2) is, 


*;(0 = A ; 
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( 4 ) 


where t DB is the time at which debonding occurs (for the particular interface), and Ay and By are parameters that 
characterize the interfacial time-dependent behavior. The typical effective interfacial behavior represented by the 
ECI model is plotted in Figure 1. Until the interfacial stress reaches the bond strength (g db ), the interfacial 
displacement is zero. Then, debonding occurs and the interfacial stress rises and then falls as the interface opens. 
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Normalized Interfacial Displacement 

Figure 1. — Comparison of the ECI model with the Needleman (Ref. 13) 
interface (Nl) model and the statistical interfacial failure (SIF) model 
developed by Robertson and Mall (Ref. 14). 

Also presented in Figure 1 is a comparison of the ECI model with two previous interfacial constitutive models 
that also admit unloading of the interfacial stress. The well-known model of Needleman (Ref. 13) (NI model) is one 
such interfacial representation. While this model does not incorporate finite bond strength, the effective interfacial 
behavior is qualitatively similar to that of the ECI model. The fact that the NI model equations are explicitly 
nonlinear in their relation between interfacial tractions and interfacial displacements, can add difficulty to the global 
solution when incorporated within micromechanics models. Further, the effective interfacial response of the 
statistical interfacial failure (SIF) model, developed by Robertson and Mall (Ref. 14), is plotted in Figure 1 as well. 
Its qualitative similarity to the ECI model and NI model is obvious. Since the SIF model, like the NI model, is 
nonlinear, Robertson and Mall (Ref. 14) incorporated a linear approximation of the SIF model (also plotted in 
Fig. 1) within their employed micromechanics model (the method of cells). 

The displacement discontinuity approach with evolving compliance interface (ECI) model available within 
NASA’s FEAMAC micromechanics analysis framework is primarily developed for modeling debonding between 
fiber and matrix materials. The FEAMAC framework and the ECI model are employed in the present investigation 
to examine crack propagation in the context of finite element simulations. Detailed comparisons of the FEAMAC 
predictions against finite element-based VCCT and cohesive element predictions under Mode-I and Mode-II are 
discussed and the advantages and disadvantages of each technique are discussed in the subsequent sections. 
Prediction of the crack path by FEAMAC is compared with XFEM predictions in the absence of experimental data. 

III. Results 

A. Simulations With a DCB Specimen 

The double cantilever beam test has widely been used for obtaining Mode-I interlaminar fracture toughness in 
composite materials. The test uses a geometry shown in Figure 2 with an initial delamination crack and the 
delamination is forced to grow by pulling apart the two plies of the specimen. To begin the comparison between the 
FEAMAC, VCCT, and the cohesive element approaches, a transversely-isotropic material and the geometry shown 
in Figure 2 is considered. The selected beam geometry has a length (L) of 4 in., width (2W) of 0.12 in. and thickness 
of 0.3 in. as used by Song et al. (Ref. 15). The beam includes two sub-laminates, each with a thickness of 0.06 in. 
and an initial crack length (C L ) of 1.15 in. between the two sub-laminates. Material properties of AS4/3501-6 
(Ref. 15) given in Table 1 are used in the plane-strain finite element models, considering unidirectional fibers 
aligned with the length direction (Y) of the specimen. 
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◄ L ► 

Figure 2. — Geometry of the DCB specimen with applied boundary conditions. 

TABLE L— MECHANICAL PROPERTIES OF AS4/3501-6 INCLUDING PARAMETERS USED WITH 
TRACTION-SEPARATION BASED COHESIVE ELEMENT SIMULATIONS (REF, 15) 


Efksi 

E 2 /ksi 

E 3 /ksi 

G 12 /ksi 

G 13 /ksi 

G 23 /ksi 

21,500 

1640 

1640 

871 

871 

522 


v 12 _ v 13 

v 23 

G IC /in. -lb/in. 2 

Gnc/in.-lb/in. 2 

Gmc/in. -lb/in. 

h 

0.30 

0.45 

0.486 

3.17 

3.17 

1.8 


a°/psi 

a°/psi 

a^/psi 

Kfpsi 

K 2 /psi 

K 3 /psi 

7800 

12600 

12600 

1.37e9 

1.37e9 

1.37e9 


The middle strip of elements shown in blue in Figure 2 in the case of the FEAMAC simulations, are represented 
with a FEAMAC micromechanics repeating unit cell (RUC) through an ABAQUS user material. Usually such an 
RUC consists of a number of subcells and is used to represent a composite material. However, for these DCB 
simulations, the material is treated as homogenized and transversely isotropic, so the RUC consists of a single 
graphite/epoxy subcell with X-Z plane of isotropy, and interfacial debonding activated at the RUC interface as 
indicated in Figure 2. The same elastic transversely-isotropic graphite/epoxy material properties are used within the 
single-subcell FEAMAC RUC. It should be noted that FEAMAC enforces periodicity conditions at the RUC 
boundaries, thus the shown RUC represents a material that is continuous in the Y-direction, even though an 
interfacial displacement discontinuity, according to the ECI model described above, is enabled in the X-direction. 
The remainder of the geometry (yellow) is meshed using bilinear incompatible-mode plane strain (CPE4I) elements 
with elastic orthotropic material model with transverse isotropy as given in Table 1. 

In the case of simulations utilizing the ABAQUS cohesive and VCCT approaches, the geometry shown in 
Figure 2 remains the same with the exception of the introduction of a surface between the two plies. For simulations 
with the cohesive approach, cohesive surface elements are employed to discretize this surface between the two plies, 
whereas simulations utilizing the VCCT approach, this surface between the two plies is defined as the possible crack 
path. In both approaches, same crack length of C L is maintained and the remainder of the model is meshed using 
bilinear incompatible-mode plane strain (CPE4I) elements with the standard ABAQUS elastic orthotropic model 
with transverse isotropy, representing the same graphite/epoxy material discussed above. 

B. Comparison of VCCT, Cohesive, and FEAMAC predictions for DCB Test 

A comparison of reaction forces predicted by the FEAMAC approach, VCCT approach, and cohesive surface 
approach is shown in Figure 3. In the case of VCCT simulation, debonding is determined entirely by the critical 
strain energy release rate Gi C reported in Table 1. For cohesive surface approach, normal debond strength g° as well 
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Figure 3. — Reaction forces predicted by VCCT, cohesive, and 
FEAMAC using adjusted G c and cr° for DCB. 
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Figure 4. — Comparison of stress component along the fiber (vertical) direction for cohesive, VCCT, 
and FEAMAC models at the beginning of crack propagation for a similar width of the geometry. 


as G IC have been used to determine the crack initiation and propagation, respectively. For the FEAMAC approach, 
this same debond strength c>° is used for the bond strength parameter as expressed in Equation (2). In addition 
to the g£ b , two other parameters critical to interfacial time-dependent behavior of the ECI model are the A t and Bj 
as reported in Equation (4). To create an interfacial debonding behavior similar to the cohesive and VCCT 
approaches, a value of 0.1 is selected for A t and 3.0 is selected for the parameter Bj. 

The stress distributions for these three models along the crack front are shown in Figure 4. As expected with the 
denser cohesive mesh, the local stress distribution near the crack- tip is relatively high compared to the VCCT and 
FEAMAC DCB models. Both VCCT and FEAMAC models use the same element size of 0.01 in. over the entire 
length of the beam. Regardless of the disagreement of local peak stresses near the crack-tip, the overall stress 
distribution shows similarity while the overall response of these three models remain consistant as evidenced from 
reaction forces shown in Figure 3. 
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TABLE 2.— TOTAL CPU TIME TAKE FOR COMPLETING 
A SIMULATION WITH FEAMAC, VCCT, 

AND COHESIVE MODELS 


Model 

Total CPU Time/s 

FEAMAC Model 

1423 

VCCT Model 

3241 

Cohesive Model 

14,806 


A 

w=0 


displacement 


Crack tip 


0.5L 


i 

L 

1 

r 

* a 0 

1" c 0 *1 


\ 

u, 


=0 


Figure 5. — Geometry of the ENF specimen with applied boundary conditions (Ref. 15). 

The time taken to complete the execution of each model is shown in Table 2. Clearly, a significant execution 
time advantage is achieved when employing the FEAMAC approach over the VCCT and cohesive approaches. The 
DCB model with cohesive elements having a maximum allowable time of 0.01 sec, took approximately ten times 
more CPU time as compared to the FEAMAC simulation with a 0.075 sec constant time increment. If a larger 
maximum allowable time increment is used for the cohesive model, the solution process becomes unstable due to 
severe nonlinearities at the verge of the crack initiation. Abaqus/Standard offers a set of stabilization mechanisms to 
handle nonlinear problems by allowing the user to specify a number of solution control parameters such as number 
of equilibrium iterations (7 0 ) and number of equilibrium iterations after which the logarithmic rate of convergence 
check begins ( I R ). Even with the current increment of 0.01 sec for the cohesive model, these convergence parameters 
had to be modified in order to stabilize the solution process. 

C. Simulations With an End-Notched Flexure Specimen 

The three-point end-notched flexure (ENF) test has been used to obtain the fracture toughness in Mode-II type 
delaminations where the initial crack is forced to grow by shearing. The dimensions and geometry used for ENF 
simulations is same as the DCB simulations, except for the applied boundary conditions (see Fig. 5), wherein one 
end is pinned and the other end is roller-supported. In contrast to the separation of laminates in the DCB test, a 
displacement at the middle of the specimen is used in the ENF simulations. 

Modeling of this geometry with FEAMAC, VCCT, and cohesive elements, is identical to the respective DCB 
simulations, that is the material data employed related to VCCT simulations (critical strain energy release rate G nc 
reported in Table 1 determines debonding) and for the cohesive surface approach, the tangential debond strength c>° 
as well as G nc are used to determine the crack initiation and propagation, respectively. For the FEAMAC model, the 
same debond strength has been used as the bond strength o t DB as expressed in Equation (3). In addition to the 
two other parameters critical to interfacial time-dependent behavior of the ECI model are the A t and Bj as 
reported in Equation (4). To create an interfacial debonding behavior similar to the cohesive and VCCT approaches, 
a value of 5x1 0” 4 is selected for A t while 0.3 is selected for the parameter Bj. 

D. Comparison of VCCT, Cohesive, and FEAMAC predictions for ENF Test 

To verify the material data presented in Table 1 related to Mode-II delamination, and identify the simulation- 
specific parameters associated with VCCT, FEAMAC, and cohesive surface approaches, ENF simulations have 
been conducted for comparison purposes. The reaction force measured at the center of the beam versus the applied 
displacement is shown in Figure 6. Clearly, all three models start debonding at approximately the same load level, 
yet the subsequent crack propagation pattern is distinctly different for the FEAMAC model as compared to both the 
VCCT and cohesive surface models. The latter two models qualitatively behave similarly with a sudden drop 
followed by a similar propagation pattern after the initiation of the crack. Note, the nature of damage progression in 
an ENF specimen is sudden and does not follow a slow crack propagation pattern as does the DCB simulation. Due 
to this sudden nature of crack propagation, it is essential to cognitive of the energies dissipated due to artificial 
damping introduced for numerical convergence purposes. Both cohesive and VCCT approaches needed significant 
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Figure 6. — Reaction forces predicted by VCCT, cohesive, 
and FEAMAC using adjusted G c and for ENF. 




Figure 7. — Initiation and propagation of a crack front between two layers. 


damping in order to maintain their numerical stability during simulations; whereas FEAMAC did not require any 
artificial damping or fine-tuning control parameters associated with the iterative convergence process. Therefore, in 
addition to execution time, the time needed to condition the solution process, i.e., experiment with numerous control 
parameters with FEAMAC is minimal, if not nonexistant, as compared to the VCCT and cohesive surface 
approches. 

E. Effect of Mesh Size on FEAMAC 

Finite element simulations of the damage and failure of composite materials have displayed a strong mesh 
dependency in crack initiation and propagation. This has been reflected on the global load-displacement response as 
well as the onset, size, and orientation of localized damage zones and crack growth in the composite materials. 
According to the implementation of FEAMAC, the stiffness of the entire element in a given direction reduces as the 
crack initiates at an integration point of an element, and progresses between two layers. Due to this stiffness 
degradation, a large crack opening (compared to a crack tip) is created and the mesh has to be fine enough to 
represent the propagating crack front (see Fig. 7). This mesh-dependence is common to other continuum damage 
mechanics approaches which rely upon stiffness degradation as their load redistribution mechanism. Consequently, 
to simulate the crack propagation as accurately as possible, the mesh must be refined along the possible crack front 
so that the debonding becomes more of a localized phenomenon, and the loss of element stiffness does not create an 
unrealistically large crack opening surface which would decrease the crack driving forces. The refined mesh used for 
FEAMAC is shown in Figure 8(a) together with the mesh used for VCCT simulations in which the mesh can be 
relatively coarser. 
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(a). FEAMAC Mesh (b). VCCT Mesh 

Figure 8. — Mesh refinement needed for FEAMAC compared to VCCT model. 


F. Uncoupled Normal and Shear in FEAMAC 

According to the ECI model in Equation (1), debonding can take place due to a normal or shear stress on any 
given face of a unit-cell. When one of these stresses exceeds the corresponding debonding strength, debonding 
initiates and the stiffness at an integration point has to degrade for both normal as well as shear. It was found that the 
current version of FEAMAC has not been implemented to account for this situation properly, in that if propogation 
of a crack were determined completely by the shear stress, the normal debond strength would still remain as shown 
in Figure 9(b), whereas conversely, if debonding initiates due to normal stress along the face of a unit-cell, it will 
still maintain a shear debond strength and the crack would continue to grow due to the normal debonding as in 
Figure 9(c). To overcome this situation a Hashin’s failure criterion has been used to couple the normal and shear 
stresses on a given face normal to the x 2 or x 3 direction as expressed in Equation (8): 
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so that both normal and shear degradation of a given interface intiate simultaneously, if the stresses at the integration 
point satisfy the condition given in Equation (8). The simulations performed after this implementation of Hashin’s 
criteria displayed the desired effect of combined normal and shear failure as the crack continues to grow in a non- 
self-similar fashion. 

G. Prediction of Crack Path with FEAMAC 

In the previous sections, DCB and ENF simulations have been used to verify the material data for simulations 
under Mode-I and Mode-II loading conditions. Additionally, these simulations have been used to identify the 
appropriate simulation-specific conditions such as viscous regularization, mesh size, control parameters related to 
convergence, etc., for VCCT and cohesive simulations, as well as the A t and Bj parameters for the FEAMAC 
simulations. In this section, the ECI model available within FEAMAC, with the newly implemented Hashin’s 
criterion for interfacial failure, is used to demonstrate the prediction of crack path, when the crack path is unknown. 

To compare the FEAMAC predictions against other available techniques such as XFEM, simulation of a three- 
point bend specimen under the plane strain conditions subjected to Mode-I loading, has been selected. The total 
length of the specimen is taken as 2 in. with a height of 0.4 in. (Ref. 17), and an initial crack length of 0.08 in. was 
employed as shown in Figure 10. 
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(b) Crack propagation due to shear debonding 



(c) Crack propagation due to normal debonding 


Figure 9. — Shear and normal debonding within a 3-point bend specimen. 



Figure 10. — Three-point bend specimen with initial crack, boundary conditions, and loading. 

An angular rotation of 0.003 radians is applied at the center of the left and the right edges of the test specimen, 
while maintaining kinematic coupling constraints between each center point and the corresponding edge. Since the 
focus of this research is on the study of applying FEAMAC for predicting crack propagation with monolothic 
materials, isotropic material properties with damage parameters reported under Table 3 are used for these 
simulations. 

Simulation of the problem with VCCT, FEAMAC, and cohesive surface approaches are very much similar to 
the simulation of DCB under Mode-I loading using each technique. Unlike VCCT and cohesive approaches, XFEM 
and FEAMAC do not require a pre-existing crack or any prior knowledge of the crack path. Since the loading is 
symmetric on the three-point bending specimen, the crack is supposed to follow the expected crack path, as shown 
in the Figure 10. In the cohesive and VCCT models, this path has to be specified as the prospective crack path and 
the models will not look for other possible crack directions. However, in the XFEM simulation, crack direction can 
change based on the specified crack propagation criteria. Any purturbation in the numerical solution process may 
cause the crack to change its direction, even though the loading, boundary conditions, geometry, and meshing are 
perfectly symmetrical (Ref. 17). Therefore, to compare the results with the VCCT and cohesive approaches, it is 
necessary to maintain the symmetry of displacement by constraining the displacements of the left and right edges 
where the angular displacement is applied. This precaution is necessary for XFEM simulations to force the 
propagation of the crack in the anticipated direction. In fact, the need to utilize constraint equations to maintain 
symmetry is a deficiency in the implementation of XFEM, since the crack should have propagated along the path 
highlighted in Figure 10 without introducing additional constraints. 
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TABLE 3.— MATERIAL PROPERTIES USED FOR THREE-POINT BEND SIMULATIONS 


E/ksi 

v 

Gic/in. -lb/in. 2 

Gnc/in.-lb/in. 2 

Gnir/in. -lb/in. 

h 

29,000 

0.3 

0.571 

0.571 

0.571 

2.284 

cr°/ksi 

a°/ksi 

a°/ksi 


25.4 

25.4 

25.4 

Only for FEAMAC 

Only for XFEM 


a s db /ksi 

A; 

Bj 

Maximum Principal Stress/ksi 

33.4 

33.4 

0.10 

0.30 

18.9 



Angle 


Figure 11. — Reaction moment predicted by FEAMAC, 
VCCT, Cohesive, and XFEM simulations. 


The reaction moment observed at the center point of the right-edge verses the applied angular displacement is 
plotted in Figure 11 for all considered techniques. Reaction values predicted at the initiation of the crack are similar 
among all the techniques, but the reaction values observed during the propagation of the crack are somewhat 
different among the models. Both VCCT and XFEM have sudden crack opening, while cohesive and FEAMAC 
predict relatively slower crack extension against the applied angular rotation. 

As expected, the crack in the VCCT and cohesive models followed the predefined crack path along the center of 
the specimen (indicated as the bonded surface in Figure 10). Since the symmetric conditions have been enforced 
during XFEM simulation, the crack continues to grow along the middle of the specimen, similar to the VCCT and 
cohesive simulations. It is worth mentioning that the addition of such constraint equations to an FEA model can 
significantly increase the computational time, especially with nonlinear material properties and large-scale structural 
problems. The displacement discontinuity approach in FEAMAC does not require additional constraints to maintain 
symmetry, and it was observed that the crack continued along the center similar to the three above techniques. Based 
on the results reported in Figure 11, it can be concluded that the ability to predict the crack path, and the resultant 
reaction moment from each technique, have been verified through the simulations with three-point bend specimen. 
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Displacement of 0.0016 in. 


Figure 12. — Eccentrically loaded three-point bend specimen with initial crack, boundary conditions, and loading. 



FEAMAC together with crack predicted 
by XFEM (in red) 


Figure 13. — Comparison of crack path predicted by XFEM and FEAMAC. 

The finite element-based VCCT or cohesive element approaches require that the path of the crack be known in 
advance, whereas FEAMAC and XFEM have the capability to predict where the crack is going to grow without any 
prior knowledge about the crack path. To demonstrate this capability with FEAMAC, an eccentrically loaded 
3-point bend specimen with the same dimensions and initial crack length as in the case described above, is used with 
displacement and boundary conditions as shown in Figure 12. 

Both the XFEM and FEAMAC approaches are used to obtain the crack propagation through finite element 
simulations. Predicted crack patterns for the given displacement are shown in Figure 13(a) for XFEM and 
Figure 13(b) for FEAMAC. Results shown in Figure 13(a) display the XFEM predicted crack path (shown in white) 
together with the partially failed elements as the crack propagates through the beam, whereas Figure 13(b) shows the 
the XFEM crack path (in red) superimposed on the FEAMAC elements that failed completely during the crack 
propagation. These two figures demonstrate a fundamental difference between the XFEM and FEAMAC 
approaches, in that the XFEM approach allows an actual crack surface to be generated by the crack front whereas 
theFEAMAC crack path is described by the elements that become near zero-stiffness elements to mimic the 
propagating crack, creating a property discontinuity. Therefore, FEAMAC will have the typical mesh sensitivity 
exhibited by continuum damage mechanics approaches where no physical discontinuites are present. 
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According to these simulation results, both approaches predict the crack path closely at the beginning of the 
crack propagation, but they deviate as the crack progresses along the specimen. Additionally, FEAMAC tends to fail 
a number of elements across the width of the crack path whereas XFEM fails only the elements that are directly 
divided by the propagating crack front. This discrepency could be associated with the debonding parameters utlized, 
the use of the Hashin failure criteria, and/or mesh sensitivity issues; all of which will be the subject of future work. 
Selection of an appropriate principal stress value for XFEM simulations has a significant effect on the prediction of 
the crack path as well as the extent of the the crack for a given displacement. Adequate mesh refinement is another 
factor to obtain a reliable prediction for crack path with XFEM. Since the entire element fails as the crack 
progresses, FEAMAC also relies on a finer mesh density to maintain a narrow crack path through the material. 

IV. Conclusions 

An extensive comparison of the displacement discontinuity approach in FEAMAC with existing approaches 
(cohesive elements, the Virtual Crack Closure Technique, and XFEM method) available in Abaqus finite element 
software has been conducted to evaluate the micromechanics-based approach and identify needed enhancements. As 
a result of these simulations, issues related to the implementation of FEAMAC as well as conceptual issues such as 
shear and normal coupling at debonding, have been identified and resolved. During this study, the FEAMAC 
approach has been used to investigate failure within materials efficiently and to explore the possibility of predicting 
non-self-similar crack paths. An eccentrically loaded 3-point bend test has been simulated with FEAMAC and the 
results are similar to XFEM predictions at the beginning of crack initiation, with some deviation from the XFEM 
predictions as the crack front progresses through the specimen under displacement control. 

Based on the experience with finite element based approaches and the micromechanics approach with 
FEAMAC, it is understood that the micromechanics based FEAMAC approach has far less sensitivity to parameters 
associated with numerical convergence. Additionally, FEAMAC offers significant advantage over the number of 
elements used compared to the cohesive approach and hence the computational time required to complete a 
simulation. This may offer an advantage over established finite element-based approaches in providing efficient and 
robust design tools for bonded composite joints. Most significantly, FEAMAC appears to have the potential to 
predict arbitrary crack paths similar to XFEM, whereas the VCCT and cohesive element approaches require 
predetermined crack paths. 

To validate fully the FEAMAC simulation results, it is necessary to compare the predictions with experimental 
data. In the absence of such data, comparison has been made to an XFEM simulation of an isotropic, eccentrically- 
loaded, three-point bend specimen. Isotropic material properties where utilized due to the difficulty in estimating a 
suitable principal stress value for XFEM, which prevented using orthotropic material properties. In contrast 
FEAMAC simulations need only the debonding strength values that have been verified through DCB and ENF 
simulations and thus easily applied to orthotropic and composite materials. 
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